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lated systems within the adiabatic-connection fluctuation-dissipation formulation of 
density-functional theory. We perform a quantitative comparison of a set of model 
exchange-correlation kernels originally derived for the homogeneous electron gas 
(HEG), including the recently-introduced renormalized adiabatic local-density ap¬ 
proximation (rALDA) and also kernels which (a) satisfy known exact limits of the 
HEG, (b) carry a frequency dependence or (c) display a divergence for small 
wavevectors. After generalizing the kernels to inhomogeneous systems through a 
reciprocal-space averaging procedure, we calculate the lattice constants and bulk 
moduli of a test set of 10 solids consisting of tetrahedrally-bonded semiconductors 
(G, Si, SiG), ionic compounds (MgO, LiGl, LiF) and metals (Al, Na, Gu, Pd). We 
also consider the atomization energy of the H 2 molecule. We compare the results 
calculated with different kernels to those obtained from the random-phase approxi¬ 
mation (RPA) and to experimental measurements. We demonstrate that the model 
kernels correct the RPA’s tendency to overestimate the magnitude of the correlation 
energy whilst maintaining a high-accuracy description of structural properties. 
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I. INTRODUCTION 


The remarkable rise of density-functional theory^ (DFT) in the last few decades owes 
much to the efficient treatment of exchange and correlation within the local-density approx¬ 
imation (LDA).^ However, known dehciencies in the LDA’s description of certain systems 
has led to the development of a hierarchy of exchange-correlation functionals of varying 
computational expense.^ At the high-complexity end of this spectrum of functionals lies the 
adiabatic connection-fluctuation dissipation formulation of DFT (ACFD-DFT),^’^ which in 
its simplest form corresponds to the random-phase approximation (RPA) correlation energy.® 
“Beyond-RPA” methods strive for an even higher level of accuracy and form an important 
and fast-developing held of research. 

ACFD-DFT provides a natural path for the improvement of the RPA via the introduc¬ 
tion of an exchange-correlation (XC) kernel f^c, an ubiquitous quantity in time-dependent 
DFT (TD-DFT).^^’^^ The homogeneous electron gas (HEG) has become a key system for 
the development and testing of new kernels through the ACFD-DFT calculation of correla¬ 
tion energies. Jellium slabs also form important test systems, through the ACFD-DFT 
calculation of their surface energies,inter-slab interaction energies,and as bench¬ 
marks for widely-used semilocal XC functionals.However the last few years have seen 
the application of the ACFD-DFT formalism to calculate the correlation energies of atoms, 
molecules and solids, with promising results. 

One such example is the “renormalized kernel approach”, which was introduced based on a 
model XC-kernel named the renormalized adiabatic LDA (rALDA).^®’^^ The rALDA exploits 
the accurate reciprocal-space description of HEG correlation provided by the adiabatic LDA 
(ALDA) in the long-wavelength limit, whilst correcting the ALDA’s unphysical behavior at 
short wavelengths. The rALDA and its generalized-gradient analogue (rAPBE) have been 
shown to yield highly accurate atomization and cohesive energies of molecules and solids. 

It is interesting to place the rALDA into the context of other HEG kernels. Many 
theoretical studies have explored the properties of the exact XC-kernel and derived certain 
limits which are not necessarily obeyed by the rALDA.Similarly, the rALDA is static, 
and apart from studies of the HEG^® there is little known about dynamical effects on ACFD- 
DFT correlation energies. Furthermore the XC-kernel of an insulator is known to behave 
qualitatively differently to a metallic system like the HEG,^® with the XC-kernel of an 
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insulator famously diverging oc l/Zc^ in the limit of small wavevectors k.^'^ In this respect it 
is important to test the validity of applying a model HEG kernel to non-metallic systems. 

This work explores the above aspects through a quantitative comparison of model HEG 
XG-kernels. Within our sample of XG-kernels we include the rALDA,^° and also a kernel 
which satishes exact limits of the HEG/^ a simple dynamical kernel/® and a kernel which 
has a divergence oc when describing an insulator.^® For each XG-kernel we use AGFD- 
DFT to calculate the correlation energy of a test set of 10 crystalline solids and evaluate the 
lattice constant and bulk modulus, which we then compare to calculations using semilocal 
functionals and the RPA, and also to experiment. We also provide a demonstrative calcu¬ 
lation of the atomization energy of the hydrogen molecule to highlight the importance of 
spin-polarization. We hnd that all of the model XG-kernels greatly improve the magnitude 
of the RPA correlation energy whilst providing a highly accurate description of structural 
properties. 

Our study is organized as follows. In section II we review AGFD-DFT and the role played 
by the XG-kernel. In particular, we describe the expected behavior of the XG-kernel for the 
HEG at certain limits (section HG), introduce our chosen set of model kernels (section HD) 
and apply them to the HEG (section HE). For inhomogeneous systems we require a scheme 
to generalize HEG kernels for a varying density; in section H G we discuss possible schemes 
and justify the choice made in this work. Section HI contains the results of our study, 
in which we discuss the calculated lattice constants and bulk moduli, absolute correlation 
energies and the H 2 molecule. Finally in section IV we summarize our results and offer our 
conclusions. 


II. THEORY 

A. Correlation energies in the ACFD-DFT framework 

Here we summarize the essential concepts of AGFD-DFT. Full derivations may be found 
in original articles"^’® or recent reviews, e.g. Refs. 6 and 39. 

In AGFD-DFT, a system of fully-interacting electrons is described by a coupling-constant 
dependent Hamiltonian H{\). The coupling constant A takes values between 0 and 1 and 
dehnes an effective interaction between electrons as Xvc, where Vc is the Goulomb interaction. 
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H{X = 1) corresponds to the exact Hamiltonian of the fully-interacting system. In addition 
to the effective Coulomb interaction H{X) contains a A-dependent single-particle potential 
Ugxt, constructed in such a way that the ground-state solution of H{X) has exactly the same 
electronic density as the ground-state solution of the fully-interacting (A = 1) Hamiltonian. 
The fixed-density path connecting A = 0 and 1 dehnes the “adiabatic connection”. Since 
H{X = 0) describes a system of non-interacting electrons with a fully-interacting density, 
is readily identihed as the Kohn-Sham potential from DFT.^’^ 

Invoking the Hellmann-Feynman theorem, integrating with respect to A along the adia¬ 
batic connection and comparing to standard DFT^ yields an expression for the exchange- 
correlation (XC) energy in terms of the operator describing density fluctuations. The 
fluctuation-dissipation theorem^® provides the link between this operator and the frequency 
integral of a response function For non-interacting electrons, = Xks, the Kohn- 
Sham response function of time-dependent DFT.^^’^^ The XC-energy is then written as the 
sum of an “exact” exchange contribution and the correlation energy Ec, where the latter 
is given by (in Hartree units): 

roo P -| 

Ec = -—J2 ds Tt Vciq)ix^iq, is) - XKsic[, is)) . (1) 

ZTT q JO JO -I 

Equation 1 has been written in a plane-wave basis so that the quantities on the right-hand 
side are matrices in the reciprocal lattice vectors G and G', and the wavevectors q belong to 
the hrst Brillouin zone. The Coulomb interaction is diagonal in a plane-wave representation 
with elements 47r/|q-|- Gp, and s is a real number corresponding to an imaginary frequency, 
u = is. 

The link between the interacting and non-interacting response functions is supplied by 
linear-response theory,which describes the behavior of density n in the presence of a small 
perturbation: 

(5n(q, u) = x^(q, cj)5ng,,t(q, uj). (2) 

The fact that yields the exact density response at all values of A allows the link to be 
made to xks through the following integral equation, 

x^(q, (^) = XKs(q, (^) + XKs(q, ^^)x^(q, ^) (3) 

where the Hartree-XC kernel f^xc has been introduced as 

fLM, w). (4) 


4 



The XC-kernel describes the change in the XC-potential upon perturbing the density, 
which is a fully nonlocal quantity in time and spaced^ 




5n(r', t') 


(5) 


By assuming approximate forms for xks and one can use equation 3 to calculate and 
thus evaluate the correlation energy with equation 1. 


B. ACFD-DFT in practice 


In a plane-wave basis, the Kohn-Sham response function has the form^^ 


XKS^'(q,*s) = 


/i^'k+q) 




IS “h '^i^k+q 


( 6 ) 


where and represent the occupation factor and energy of the Kohn-Sham state d’z/k, 
while the pair densities nyk,i^'k+q(G) are matrix elements of plane waves, (t/’iyk|e“**'‘^’'“‘^^'“'|d’i/'k+q)- 
n is the volume of the primitive unit cell, and the factor of 2 assumes a spin-degenerate 
system. From equations 1 and 6 the benehts of the ACFD-DFT are not very obvious; to 
construct yxs we require ip, which means solving the Kohn-Sham equations and thus already 
obtaining the correlation energy. Furthermore, to solve the integral equation (3) we require 
the XC-kernel jPc, which arguably is even more complicated than the XC-potential v^c due 
to its frequency dependence. 

However the attraction of ACFD-DFT is that even setting jPc = 0 yields both a nonlocal 
description of exchange and a nontrivial expression for the correlation energy, namely that 
obtained from the RPAd 


1 _ roo 

I ds Tr [In {1 - nc(q)yKs(q, +'yc(q)yKs(q, *«)] • (7) 

The RPA has been applied across a wide range of physical systems^^^^® and found to give a 
markedly improved description of nonlocal correlation effects. Equation 7 is usually applied 
as a post-processing step to a DFT calculation, analogous to GohFo corrections to band 
gaps.^°’^^ 

Based on the success of the RPA it may be hoped that the description of correlation 
might be further improved by using more sophisticated approximations for fxc- While it 
turns out that the adiabatic local-density approximation (ALDA) offers little improvement,^^ 
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nonlocal, dynamical and/or energy-optimized approximations for f^c have been found to 
correct deficiencies of the RPA when calculating the correlation energy of the homogeneous 
electron gas 

We note that a non-self-consistent application of the ACFD-DFT formula might suffer 
from a dependence on Vxc, the exchange-correlation potential used to construct the orbitals 
forming xks- In this respect, a self-consistent scheme is attractive and the subject of current 
researchP^^®^ However here we do not include any self-consistency and treat f^c as a quantity 
to be optimized independent of the Vxc used to generate Xks-^® 

C. XC-kernels from the homogeneous electron gas 

In the same way that the HFG is used to generate approximate XG-potentials, it also 
forms a natural starting point for approximate XG-kernels. Here we review some properties 
of the exact XG-kernel of the HFG. 

1. Definitions 

The analogue of equation 3 for the HFG is: 

X^{k,uj) =Xo{k,uj)+ Xo{k,uj)fHxc’^{k,uj)x^{k,uj). (8) 

All quantities appearing in this equation are scalars, with /c = |G -|- q|. Xo is the Lindhard 
response function, with occupation numbers equal to 1 for plane-wave states below the 
Fermi level and zero otherwise.®^ As demonstrated in the appendix of Ref. 58, in the case that 
equation 3 is applied to the HFG the Lindhard and Kohn-Sham response functions coincide, 
so that the quantity appearing in equation 8 must also match its counterpart 

in equation 3. Therefore the HFG forms a rigorous test ground for model XG kernels. 

For simplicity we denote as The local held factor G{k,u) is related to 

/™^(A;,a;) as /™^(A;,a;) = —Vc{k) G{k,u), so that equation 8 dehnes G in terms of the 
response functions x^^^ A potential source of confusion is that another held factor 

Gj can be found in the literature which has a diherent dehning equation.For Gj, 
the Lindhard response function appearing in equation 8 is modihed such that the occupation 
numbers of each plane-wave state are calculated using the many-body wavefunction of the 
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fully-interacting system.^® Therefore is also a distinct quantity. However the kernels 

investigated in this work were derived based on the equivalence of equations 3 and 8 for the 
HEGZ® so it is G which is of interest in the current work. 


2. Exact limits 


The local-held factor G (and thns have been the subject of many theoretical studies 

(c.f. section IIIC of Ref. 33), and their behavior at certain limits is known exactly. First, in 
the long wavelength {k —)■ 0) and static (ca = 0) limit, the HEG XG-kernel rednces to the 
adiabatic local-density approximation (ALDA): 

47rA 




ALDA _ 
xc 


kp 


where 


^ _ 1 _ kl (f{nec 
4 47r dri^ 


( 9 ) 


( 10 ) 


kp = (37r^n)^/^ is the Fermi wavevector for the HEG of density n, and Ec is the correlation 
energy per electron. The two terms in equation 10 correspond to the exchange and correla¬ 
tion contribntions to the ALDA kernel. Eqnation 9 can be seen either as a consequence of 
the compressibility snm rnle,^^ or more simply by noting that the ALDA should be exact in 
describing the HEG response to a nniform, static held.^® 

Remaining in the static case, but considering small wavelengths {k cx)) yields^^’®^ 

47rR 47rG 


fxc (/c ^ cx,cu = 0) = 


k^ 


kp 


( 11 ) 


whilst in the long wavelength, high frequency limit {k = 0,u ^ oo) is given by^^ 

47rD 


= 0,u^oo) = 


kl 


( 12 ) 


Although we have not written it explicitly. A, B, G and D depend on the density of the 
HEG, or eqnivalently on the Fermi wavevector or Wigner radius = (3/47rn)^/^. Practically 
A, G and D can be obtained from a parameterization of the HEG correlation energy £c, while 
B additionally reqnires the momentnm distribntion and on-top pair-distribntion fnnction of 
the HEG.®^ In this work we use the parameterization of Sc and B from Refs. 63 and 64 
respectively. 
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3. Intermediate k values 


In addition to these limits, calculating the correlation energy from equation 1 requires an 
expression for across all k and u. Ref. 64 provided important insight into the HEG 

XC-kernel with diffusion Monte Carlo calculations of ca = 0) for a range of k vectors 

and densities. Interestingly, one of the conclusions of the work was that in this static limit 
and for k < ^kp, could be well-approximated by taking its /c = 0 limiting form (i.e. 

the ALDA, equation 9). 

4 . The large-k limit and the pair-distribution function 

According to equation 11, the /c —)■ cx),a; = 0 limit of is a constant, and thus 

G{k,u = 0) diverges as k"^. This property appears to have alarming consequences for the 
pair-distribution function: in Ref. 65 it is shown that a G{k) which diverges for large k will 
yield a singular pair-distribution function at the origin (see also Refs. 27,66). The crucial 
point to note however is that the relationship in Ref. 65 assumes a frequency-independent 
G{k), which is not the same as the frequency-dependent G{k,u) evaluated af a; = 0, which 
equation 11 describes. This point is discussed further in Ref. 58. 

Nonetheless, a frequency-independent model kernel satisfying the exact HEG limits for 
small and large /c at ca = 0 (e.g. the CDOP kerneP^) will have a badly-behaved pair- 
distribution function, as demonstrated below and in Fig. 2 of Ref. 29. In this respect, and 
given that the correlation energy is calculated as an integral over frequency, a frequency- 
averaged G{k) may be a better starting point for model kernels than G{k,u = 0). Below 
we compare kernels which do or do not satisfy equation 11. 

D. Model XC kernels 

Having introduced the relevant quantities for the HEG, we now describe the XG-kernels 
investigated in the current work. The XG-kernels are plotted in reciprocal or real space 
in Figs. 1(a) and (b). The reciprocal and real space XG-kernels are related by the Fourier 
transform 


VcW, At, w) ^ J dRe R, W) 


(13) 






FIG. 1. Comparison of model HEG XC-kernels, plotted in (a) reciprocal and (b) real space for 
rg = 2. In (b) the XC kernels are divided by the Coulomb interaction Vc, and the inset provides 
a zoomed image close to fxc = 0. Apart from the JGMs all the XC-kernels are short-range, while 
(apart from CDOPs) at small R the XC-kernels cancel the Coulomb interaction such that fuxc 
vanishes. We have omitted the CDOP kernel in (b), since it matches CDOPs apart from a 5- 
function at R = 0 (equation 18). The CPd kernel was evaluated at an energy of 2 Hartrees and 
the JGMs kernel at a band gap of 3.4 eV. In (c) we plot the wavevector-resolved correlation energy 
(equation 29) at = 4 compared to the exact^® result obtained from the parameterization of the 
correlation hole given in Refs. 67 and 68. In (d) we plot the difference in calculated correlation 
energies with a parameterization®^ of Monte Carlo calculations®® of Ec- 

with R = |r — r'l. 


1. The rALDA kernel 

The renormalized adiabatic local density approximation (rALDA)®® XC-kernel is given 
by 


k) 

where the Heaviside function 6{x) 
is chosen as 


e{kc - k)^ + e{k - ke)^ 


(14) 


1 for X > 0 and zero otherwise. The cutoff wavevector 


kc = kp/^/A 


(15) 
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In previous applications of this kernel^^’^^ the coefficient A defined by equation 10 was 
replaced by 1/4, corresponding to omitting the correlation contribution. In this work we 
shall refer to this exchange-only kernel as the rALDA kernel The special label 

(rALDAc) refers to the kernel calculated including both the exchange and correlation 
contributions in equation 10. We note that the rALDAc kernel coincides with that of Ref. 61 
with B = 1 and C = 0. 

obeys the exact k 0, lj = 0 limit for the HEG (equation 9). Furthermore, 
both and mimic the HEG kernel®^ in displaying small variation for wavevec- 

tors below ^ 2kp- At larger wavevectors both kernels correspond exactly to the Goulomb 
interaction with opposite sign, such that the corresponding Hartree-XG kernels Jhxc vanish 
for k > kc- In real space the kernel has the form 

frALDA/ O'! _ ^ 

Jxc 

with the Fourier transform of the Heaviside functions leading to decaying oscillations 
[Fig. 1(b)]. At small R the XG-kernel diverges as —1/R, yielding a Hartree-XG kernel 
which is hnite at the origin.^*’ 



kcR 


sin X [sin ( kcR) — kcR cos ( kcR)] 


-dx 


X 


{k,Rf 


( 16 ) 


2. The CD OP kernel 


The kernel introduced by Gorradini, del Sole, Onida and Palummo (GDOP) in Ref. 17 
has the form 




dTTQ; 


A 1 


+ 


47rG 


(17) 


kp J kj, [g + {k/kpy] kip 

where g = B/{A — C), and a and jd are density-dependent fitting parameters chosen to best 
reproduce the local held factor G{k,u = 0) obtained from the QMG calculations of Ref. 64. 

Uniquely among the kernels considered here, the GDOP kernel obeys both the /c —)■ 0 
and k ^ oo limits of the HEG at ca = 0, equations 9 and 11. However as noted above 
the short-wavelength C term causes the pair-distribution function to diverge.^® In Ref. 29 
a simplihed kernel which avoids this divergence was obtained from equation 17 by setting 
G = 0. We shall also investigate this kernel in this work, labeled GDOPs (/A^^^)- 
The real-space form of the GDOP kernel is: 




/tt W 


V/9 / 

[ 2/3 J 




( 18 ) 
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We note that the C term in equation 17 produces a 5-function in real-space, while for small 
R (excluding the 5-function) the XC-kernel diverges as —B/R, such that the Hartree-XC 
kernel is still divergent as (1 — B)/R. 


3. The CP kernel 


A kernel with a simpler functional form was introduced by Constantin and Pitarke (CP) 
in Ref. 16: 


fxc (n, k) = - 


dTT 


(1 - ) 


(19) 


Here kq = A/kp, which ensures that the HEG /c —)■ 0, ca = 0 limit is satisfied. Like the 
rALDA kernels, at large wavevectors {k) cancels the Coulomb interaction so that Jhxc 
vanishes. The CP kernel possesses a compact form in real space in terms of the error 
function: 


fxc {n,R) = -- 


erf 


R 




( 20 ) 


As R —)■ 0, diverges as —1/R and thus yields a hnite Hartree-XC kernel in this limit. 


4- The CPd dynamical kernel 


With its simple form, the CP kernel is an ideal starting point to explore more complex as¬ 
pects of fxc, such as its frequency dependence. In Ref. 16, a dynamical kernel was introduced 
(CPd) by replacing kq appearing in equation 19 with i.e. 


^CPd 
J XC 


(n, k, u) 


'dvr 


(1 - 




) 


( 21 ) 


where for imaginary frequency u = is, 





-|- as 4- cs^ 
l + s2 


( 22 ) 


In Ref. 16 the coefficient c = D/A was chosen to correctly reproduce the k ^ 0, u ^ oo 
limit of the HEG (equation 12), while the relation a = 6i/c was found to give a good £t 
to the correlation energy calculated for the HEG using We note that the CPd kernel 

varies non-monotonically with frequency in the /c —)■ 0 limit. 
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5. The JGMs kernel 


The limits of the exact kernel discussed in section IIC 2 were derived for the HEG, which 
is metallic. However, the XC-kernel of a periodic insulator is known to behave differently, 
diverging oc l/k"^ in the /c —)■ 0 limit.This limit has been found to play an essential role 
in the TD-DFT calculation of excitonic effects in optical spectra, leading to the development 
of kernels which exhibit the same l/k‘^ divergence.Here we focus on the “jellium with 
gap” model (JGM) kernel of Ref. 38, which has the useful property of reducing to a model 
HEG kerneh® when applied to a metallic system. 

The JGM kernel of Ref. 38 was derived based on two steps. First, the theoretical argu¬ 
ments of Ref. 72 were used to connect the /c —)■ 0 limit of f^c to e, the dielectric function, as 
fxc{k —)■ 0) = —47r/[A;^(e — 1)] (a similar relation was found empirically in Ref. 70). Then, 
the model dielectric function of Ref. 73 was used to relate e to the band gap of the material 
Eg, as e — 1 = Ann/Eg. The same power dependence may be found for other model dielectric 
functions, e.g. the Penn model,and essentially follows from the /-sum rule. Gombining 
these relations places a requirement on the model kernel that it diverges as —a/k'^ in the 
small k limit, where a ^ Eg/n. 

In Ref. 38 the JGM kernel was constructed to satisfy this divergence, based on a modihed 
GP kernel introduced in Ref. 16. Following their approach, starting from the unmodihed 
GP kernel we introduce a simplihed JGM (JGMs) kernel for a system with a gap 

as 


f^?^%n,k,Eg) 


_ n _ g-K0fc2g-S2/(47rn)N 

. k'^ 


(23) 


has many of the properties of the JGM kernel introduced in Ref. 38. For systems 
with a band gap, the XG-kernel diverges as —a/k"^ at small k. For Eg = 0, reduces to 

fxc (equation 19), while for Eg ^ oo —)■ —Vc, yielding a vanishing correlation energy. 

Indeed the JGMs kernel differs only from the JGM kernel in its behavior at large k, with 
the JGM kernel correctly reproducing the ca = 0 limit of the HEG for Eg = 0 (equation 11); 
concomitantly the JGM kernel has a diverging pair-distribution function. By introducing 
the JGMs kernel we can study the effect of the —divergence without any additional 
complications potentially arising from a badly-behaved pair-distribution function. 
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The real-space form of the JGMs kernel is 

= 1 - erf j . (24) 

Equation 24 and Fig. 1(b) emphasize a unique property of the JGMs kernel: it is long range. 
As a result, at large R the Hartree-XG kernel does not reduce to the bare Goulomb kernel 
but rather to an interaction weakened by a factor exp[—Eg/(47m)]. 

6 . Coupling-constant dependence 

Evaluating the integral over A in equation 1 requires the f^c kernel at an arbitrary coupling 
strength. We use the analysis of Ref. 15 to link to the fully-interacting kernel through 
the relation 

= A"^ f^c{n/X^,k/X,uj/X‘^). (25) 

The scaling of the density can be equivalently stated as Ar* or kp/X. The (exchange-only) 
rALDA XG-kernel has the useful property that 

For the JGMs kernel, we have an additional parameter Eg. For simplicity, we employ a 
scaling 

k, Eg) = X-^ /i,«^^(n/A^ k/X, Eg/X^-% (26) 

equivalent to treating Eg/n independent of A. 


7. Analogy with range-separated RPA 


It is interesting to draw comparisons with RPA methods based on the concept of range- 
separation. First we trivially relabel the Hartree-XG kernel as an effective interaction 
Ueff, i.e. Ueff = Vc + fxc, uotiug from Fig. 1(b) that for most of the kernels Ues goes to zero 
at R —)■ 0 and tends to the full Goulomb interaction at large R. Now, specializing to a 
static model XG-kernel which scales linearly with coupling constant, = -^/xc(q), we 

can partition the correlation energy in equation 1 into two contributions Ec = 




ds Tr[ln{l - nefr(q)XKs(qG-s)} +neff(q)XKs(qG'S)^ 


F. = 


rl roo P -| 

dX ds Tr /,,e(q)(x^(q,fs) - XKs(qGs)) 

„ JO JO ^ 
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Equation 27 matches the RPA expression for the correlation energy (equation 7) with Vc 
replaced with Ueff- Similarly, equation 28 matches the expression for the full correlation 
energy (equation 1) with Vc replaced by —fxc, which is generally a short-ranged interaction. 

We can focus further on the specihc example of the CP XC-kernel (equation 19), noting 
that this kernel can be linearized in A by neglecting the correlation contribution to A in 
equation 10. Then, Ues = eTf{fiR)/R, with the “range-separation parameter” p determined 
by the density through fi = kp ^ l.Q/r^. This effective interaction is often found in the 
range-separated RPA^® with p of order unity. 

We stress that equations 27 and 28 are exact for any kernel which obeys 
Since most of the XC-kernels under study here do not obey this relation, we have not 
explored equations 27 and 28 further in this work. However for XC-kernels linear in A 
(e.g. the rALDA) there may be a computational advantage in calculating and 
separately. requires only xks and not so for a given kernel it can be obtained at 
exactly the same cost as the RPA correlation energy. In fact since the effective interaction 
Ueff generally vanishes at large wavevectors, can be expected to avoid the basis-set 
convergence problems of the RPA recently highlighted in Ref. 78. Meanwhile it may be 
possible to exploit the short-range character of f^c to reduce the computational cost of 
calculating from equation 28. 

E. Calculating HEG correlation energies 

A standard test of model HEG kernels is to calculate the correlation energy per electron 
Ec from equations 1 and 8. For a given density r^, Ec can be resolved as an integral over k 
as^^ 

roo 

ec= £c{k)d{k/2kF). (29) 

Jo 

The quantity £c{k) can be compared to the Fourier transform of a suitable parameterization 
of the “exact” correlation hole obtained from Monte Carlo calculations.^®’®^’®® Alternatively 
one can compute Ec over a range of densities and compare to the parameterized result.®® 
These comparisons are made in Figs. 1(c) and(d) respectively. 

The HEG analysis has been performed a number of times®®’^®’®® so we only summarize the 
key points. The RPA correlation is too negative, while including any of the XC-kernels brings 
Ec to within 0.1 eV of the exact result. Considering the wavevector decomposition in Fig. 1(c) 
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we find that the dynamical CPd kernel provides the best description of the correlation hole 
at this density (r^ = 4), but below ~ l.bkp there is very little difference between any of 
the XC-kernels and the exact result. Indeed the ALDA (not shown) also provides a good 
description of the correlation hole at these wavevectors. At larger wavevectors, differences 
begin to emerge between the kernels, with the CP kernel becoming too negative, the CPd 
and CDOPs kernels closely following the exact result, and the other kernels too positive. 
The rALDA kernels are abruptly cut off at edk) = 0, while the CDOP kernel acquires a 
slowly-decaying positive contribution. The latter behavior is observed to a greater extent in 
the ALDA, and originates from the locality of the kernels. 

Over the full range of densities [Fig. 1(d)], we hnd that the calculated correlation energy is 
slightly too positive with the CDOP kernel and too negative with the CP and CDOPs kernels. 
Interestingly, CDOPs is closer to the exact result than CDOP, illustrating that removing 
the part of the CDOP kernel which causes the pair-distribution function to diverge^® slightly 
improves the correlation energy. The rALDA kernels fall closest to the exact result across a 
wide range of densities, and the CPd kernel also provides a good description of the correlation 
energy. Comparing rALDA and rALDAc, we see that removing the correlation contribution 
from A in equation 10 decreases the correlation energy per electron by less than ~0.02 eV 
across a range of densities. 


F. Coupling-constant averaged pair-distribution function 

Clearly all of the considered kernels greatly improve the correlation energy of the HEG 
compared to the RPA. The common characteristics shared by the kernels are that they satisfy 
the exact /c —)■ 0, ca = 0 limit of the HEG (except the rALDA, which neglects the correlation 
contribution in equation 10), and that they decay for wavevectors above 2kF- This decay 
is essential to an accurate description of the energetics of the HEG, with the ALDA (which 
does not decay at large k), yielding a correlation energy which is too positive.However 
the fact that we only observe small variations between the kernels considered in Fig. 1(d) 
indicates that the precise form of this decay is less important. 

It is however interesting to consider the coupling-constant averaged pair-distribution func¬ 
tion ^c(-R),®^ obtained as the Fourier transform of ec{k) multiplied by n/[Inkp]-^^ The pair- 
distribution function gdR) is obtained from the derivative of gdR) with respect to r^, and 
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FIG. 2. Coupling-constant averaged pair-distribution function gc{R) calculated at = 2 for the 
model HEG XC-kernels (see Fig. 1 for color code), compared to the “exact” parameterization given 
in equation (36) of Ref. 67 (see also Ref. 68). The inset shows the zoomed region around gdR = 0). 

the exact gc and gc both satisfy cusp conditions such that their slopes at i? = 0 are gener¬ 
ally nonzero.®®’®^ Concentrating on gdR), we note that in order to describe a cusp in real 
space we require Fourier components [i.e. nonzero sdk)] at large k. Indeed the analysis of 
Ref. 65 hnds that a frequency-independent kernel must decay as —■j/k‘^, where 7 < dvr, i.e. 
the Hartree-XC kernel retains a 1/fc^ term at large k. By considering Fig. 1(a) we see that 
the rALDA, CP and CPd kernels all decay as —such that their Hartree-XC kernels 
vanish, so that edk) also quickly tends to zero at large wavevectors [Fig. 1(c)]. Thus these 
kernels cannot describe the cusp. 

To illustrate this behavior, in Fig. 2 we plot gdR) calculated at = 2 for the different 
kernels, compared to the RPA and to the parameterization of Refs. 67,68. It is clear that the 
coupling-constant averaged pair-distribution functions calculated for the rALDA, CP and 
CPd kernels are far softer than those calculated for the RPA and CDOPs, whose Hartree- 
XC kernels decay oc l//c^. Meanwhile as noted above the local C term of the CDOP kernel 
causes a divergence in gdR)- 

The slope of coupling-constant averaged pair-distribution function calculated in the RPA 
is too steep, while it is improved for CDOPs. In Ref. 29 it was also found that gdR) 
calculated for CDOPs agreed well with the exact result. The good performance of the 
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CDOPs kernel for calculations of the pair-distribution function might have been anticipated 
from the fact that the coefficient B appearing in equation 17 itself is determined from 
g{R = 0).®"^ The above analysis illustrates how the precise large-/c behavior of a kernel 
affects its description of the cusp of gdR), despite playing a lesser role in the calculation of 
energetics. 

G. Applying HEG kernels to inhomogeneous systems 

In order to calculate the correlation energy of an inhomogeneous system through equa¬ 
tion 1 we require f^c evaluated in a plane-wave basis, which in general is constructed from 
the real-space kernel through 

i / dr [ (30) 

V Jv Jv 

where V is the volume of the entire crystal, consisting of Ng replicas of the unit cell of volume 
n. The question is how to incorporate into this formalism a model (m) kernel which has the 
form /™(n, k,u) or /™(n, |r —r'|,a;). In the case that the system is homogeneous (n(r) = n), 
we simply make the substitution fxc{r,r',u) f^{n, |r — r'|,a;) to get a diagonal kernel, 

gr = Saa' fZin, Iq + GU). (31) 

Alternatively, if the model kernel is fully local (independent of k, e.g. the ALDA), it is 

natural to choose the local density to construct the kernel, and obtain 

However for nonlocal kernels and inhomogeneous systems it is not obvious how one should 
construct /^.^(r, r', ca), except for two requirements. First, an arbitrary model kernel should 
be symmetric in r and r'T® 

f^c{r,r',u) = f^c{r',r,u). (33) 

Second, for the JGMs kernel, we require that in the g —)■ 0 limit the head of fxc (i-e. 

G = G' = 0) diverges as 1/g^ while the wings (G 7 ^ G' = 0) diverge no faster than 1 /g 

(Ref. 37). As shown below, this second requirement turns out to exclude previous schemes 
used in ACFD-DFT calculations, which we now briefly review. 
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1. Density symmetrization 


A symmetric kernel can be obtained by making the following substitution into equa¬ 
tion 30: 

f^c{r,r',uj) ^/™(5[n],|r-r'|,a;). (34) 

Here, iS is a functional of the density symmetric in r and r', whose possible forms span a wide 
range of complexity.Refs. 22 and 31 take a two-point average, iS[n] = l/2[n(r) -|-n(r')], 
which has an intuitive interpretation. 

A disadvantage of the symmetrization in equation 34 is that in order to evaluate equa¬ 
tion 30, it is necessary to work with a real-space representation of the kernel and perform an 
integration over the entire volume of the crystal. For short-range kernels this integral can be 
converged by sampling over a number of unit cells,but for a long-range kernel like 
[decaying as —a/{47iR)] the required sampling might be prohibitively large. Also, since the 
kernel is constructed in real space, it is not obvious how to control the 1/q divergences of 
the JGMs kernel in reciprocal space. Finally we note that the 1/R real-space divergence of 
the kernels leads to slow convergence with the real-space grid used to evaluate the integral 
in equation 30.^^ 

2. Kernel symmetrization 

An alternative approach followed in Ref. 29 is to start from a nonsymmetric form of fxc, 
which we label 

/“"(r.r'.o,) =/,^(n(r),|r-r'|,a,). (35) 

Inserting into equation 30 gives 

i |q+ G'|,c), (36) 

A symmetric kernel can then be obtained by averaging with its Hermitian conjugate, 

i.e. 

i . (37) 

This procedure can be seen equivalently^® as inserting the symmetric combination 1/2 [/^®(r, r', ca)-!- 
,r,u)] into equation 30, and therefore corresponds to a two-point average of the ker¬ 
nel; in the case that the kernel is linear in density, this scheme is equivalent to averaging 
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the density. Equation 36 has the computational advantages that the integral is over a 
single unit cell and requires only that the density (and not the kernel) is represented on 
the real-space grid. However, considering the JGMs XC-kernel, upon inserting equation 23 
into equation 36 and performing the average of equation 37 we are left with a matrix whose 
wings diverge oc 1/g^ as g —)■ 0, not 1/g as required. Therefore equation 37 is unsuitable for 
the current study. 

3. Wavevector symmetrization 

In order to correctly deal with the JGMs kernel while retaining some of the computa¬ 
tional advantages of equation 36, we follow the approach of Ref. 38 and symmetrize the 
wavevector appearing in the right hand side of equation 36 with the substitution |q-|- G'| 

I q J- G11 q -|- G^ I; 

(q, ^ (^n(r), ^|q + G| |q + G'|, . (38) 

Like equation 36, equation 38 requires the integral over the unit cell only and deals with 
the reciprocal-space form of the kernels. Using equation 38 to construct the JGMs ker¬ 
nel yields a matrix whose head and wings diverge in the g ^ 0 limit as 1/g^ or 1/g, as 
required, while the hermiticity of /^‘^'(q, ca) automatically satishes the symmetry require¬ 
ment (equation 33). Trivially, the averaging scheme will reproduce the equations 31 and 32 
when applied to systems with a homogeneous density or a local kernel, and furthermore the 
diagonal (G = G') elements coincide with those calculated with the two-point average of 
the kernel, equation 37. 

Since this work is concerned with the comparison of a large number of kernels, we have 
elected to use equation 38 on the grounds that it is relatively efficient, and can deal with 
the divergences in the JGMs kernel correctly. However the physical interpretation of the off- 
diagonal elements arising from the wavevector-symmetrization is not transparent. Although 
a two-point scheme also suffers from limitations (e.g. the two-point kernel has no knowledge 
of the medium lying between r and r'), it still remains a more intuitive procedure. The fact 
that we have to invoke an averaging scheme at all is an undesirable consequence of using 
HEG kernels to describe inhomogeneous systems. In reality the use of different schemes can 
only be justihed through testing and comparison with experiments or other calculations. 
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such as that performed in Refs. 29-32 and here. 


H. Computational details 

All calculations in this work were performed using the GPAW code.®^ The Kohn-Sham 
states and energies used to construct the response function (eqnation 6) were calculated 
using the local-density approximation to DFT^’^’®^ within the projector-augmented wave 
(PAW) framework.®^ We used 6x6x6 and 12 x 12 x 12 unshifted Monkhorst-Pack®^ 
meshes to sample the Brillouin zone for insulators and metals respectively, and constructed 
the occupation factors for each Kohn-Sham state using a Fermi-Dirac distribntion fnnction 
of width 0.01 eV. For H, He and H 2 we used a simulation cell of 6x6x7 and F-point 
sampling. 

When calculating Ec the wavefnnctions were expanded in a plane-wave basis set np to 
a maximum kinetic energy of 600 eV. Following previous studies,^^’^® we used the frozen- 
core approximation bnt included semicore states for some elements.®^ We note that norm- 
conservation was not enforced in the generation of onr PAW potentials, while it is reported 
in Ref. 78 that including norm-conservation might increase the magnitude of the RPA corre¬ 
lation energy and decrease the calculated lattice constants for certain materials. As a resnlt 
some care should be taken in making comparisons to experiment, although we expect calcu¬ 
lations including the XC-kernel to be less sensitive to this convergence issue (section IID 7). 

For the matrices representing the response function and kernel, we used a lower plane- 
wave cntoff Ecut of 400 eV (300 eV for Na and H 2 ), and a q grid matching the Brillonin zone 
sampling of the ground-state calculation. We truncated the sum-over-states appearing in 
eqnation 6 at a nnmber of bands eqnal to the nnmber of plane waves describing the response 
function, e.g. ~700 for Si. Within this approximation, the following extrapolation scheme 
is commonly used for the RPA correlation energyd®’^® 

Ef^^iE,^,) ~ ^ 00 ) + KE;^/\ (39) 

In Ref. 31 it was proposed that the same expression can be applied to the correlation 
energy calculated with the rALDA kernel. We have tested this expression for each of the 
kernels in section IID for a set of 10 materials (see section III A). As an example, in Fig. 3 
we plot the correlation energy per electron calculated for MgO as a function of E~^J‘^. As 
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FIG. 3. The correlation energy evaluated per electron for MgO at a lattice constant of 4.23 A, 
using different approximations for the XC-kernel (see Fig. 1 for the color code). Scut is the plane- 
wave cutoff used for the matrices representing the response function and XC-kernel. The circles 
represent calculated data points, and the lines are fits from equation 39. 

demonstrated by the straight lines, equation 39 apparently gives a good description of the 
correlation energy calculated for Scut >200 eV for all of the kernels. We have observed 
the same behavior across the combinations of materials and kernels. Therefore in order to 
facilitate comparison across the entire test set we will apply equation 39 for all XC-kernels. 
We point out that the correlation energy tends to converge faster (shallower lines in Fig. 3) 
when a nonzero f^c is used, and for calculating structural properties the extrapolation is 
often unnecessary. 

Constructing the XC-kernel with equation 38 is not straightforward due to the dual de¬ 
pendence on G and G'. Our current implementation distributes the rows of the f^c matrix 
among processors before evaluating the integral in equation 38. In the future it may be ap¬ 
propriate to improve performance through an interpolation scheme, as for Refs. 29 and 86. 
On the other hand, for the systems studied here the time taken to construct the kernel is 
small compared to that spent constructing the response function xks and inverting equa¬ 
tion 3. When constructing the kernel, we use the PAW all-electron density to be consistent 
with previous work.^^ The 1/q^ divergence of the Coulomb interaction and JGMs kernel was 
treated within the scheme described in Ref. 41. 
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Having obtained we evaluated equation 1 by numerical (Gauss-Legendre) integration 
over the coupling constant {N\ = 8 points) and frequency (iV^^ = 16 points, using a logarith¬ 
mic mesh).^® By virtue of the scaling relation (equation 25) we must construct the rALDA 
kernel once, a general static kernel Nx times and a dynamical kernel NxN^^ times; hence 
there is a prefactor of ~ 100 applied to computing compared to 

To calculate structural properties, we evaluated the total energy + E^ + Ec 

for seven lattice constants centered around the experimental value and £t the values to the 
Birch-Murnaghan equation of state.We used higher plane-wave cutoffs of 800 eV (900 eV 
for MgO, LiCl and LiF) to evaluate the LDA energies and E^, and used a denser sampling of 
the Brillouin zone combined with the Wigner-Seitz truncation scheme described in Ref. 88 
to calculate E^. We typically obtain converged exchange energies for insulators with a 
sampling of 10 x 10 x 10, while metals require a denser sampling®® (e.g. 20 x 20 x 20). 
Since the bulk modulus is constructed from derivatives of the energy, it is rather prone to 
numerical error, to the extent that different code implementations of the same method can 
yield different results.®® In this respect one should attach more significance to the calculated 
lattice constants than bulk moduli, since the former are more robust quantities. However 
even for the bulk moduli one expects a reduction in error when comparing different XC- 
kernels within the same computational framework. 


III. RESULTS AND DISCUSSION 

A. Lattice constants and bulk moduli 

We have selected a test set of 10 materials, consisting of 3 tetrahedrally-bonded semi¬ 
conductors (diamond C, Si and SiC), 3 ionic compounds (MgO, LiCl and LiF) and 4 metals 
(Al, Na, Cu and Pd). For each material, we used the XC-kernels introduced in section HD 
to calculate the lattice constant and bulk modulus. Here we compare these results to those 
obtained from DFT (in the LDA or from the generalized-gradient PBE functional®®), the 
RPA, and to the experimental values tabulated in Ref. 46. 
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FIG. 4. Lattice constants and bulk moduli calculated with different XC-kernels for C, MgO and 

— 3/2 

Al, vs ) where F^cut is the plane-wave cutoff of the response function and XC-kernel matrices. 

The green and blue horizontal lines give the values calculated with the LDA and obtained from 
experiment, respectively, where the experimental data were tabulated in Ref. 46 (note the LDA 
lattice constant/bulk modulus for Al [3.99 A/83 GPa] is off the scale). The values at infinite cutoff 

_o /o 

= 0) were calculated from the correlation energies extrapolated from equation 39. 

1. General trends 

Figure 4 shows the lattice constants and bulk moduli calculated for C, MgO and Al as a 
function of a quantity inversely proportional to the number of plane-waves describing 

the response function xks and XC-kernel f^c (c.f. Fig. 3). The quantities at = 0 were 

calculated from Ec extrapolated to inhnite F^cut using equation 39. We also show the values 
obtained from the LDA and experiment as horizontal lines. 

There are three key observations to be made from Fig. 4. First, for non-metallic systems 
the rALDA, rALDAc, CDOPs and CP kernels yield almost identical results, which in turn are 
very similar to the RPA. Second, the JGMs, CPd and the CDOP kernels (which respectively 
are long-range, dynamical or have a local term) display distinct behavior. For instance the 
JGMs kernel predicts smaller lattice constants and larger bulk moduli than the other XC- 
kernels. Finally, all of the XC-kernels show faster convergence with respect to E^ut compared 
to the RPA, as found for the correlation energy (Fig. 3). 

Keeping the above points in mind, we extend this analysis to the full test set and consider 
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FIG. 5. Percentage deviation from experiment of calculated lattice constants and bulk moduli for 
the test set of 10 materials. The values used to construct the plots are presented in Table I. Each 
line corresponds to a different approximation for fxc- 


each kernel in turn. The entire dataset is given in Fig. 5 and Table I. 


2. LDA, PBE and RPA 

The LDA typically underestimates lattice constants and overestimates bulk moduli, while 
PBE displays opposite behavior. For tetrahedral semiconductors the LDA is difficult to beat, 
and is by far the most computationally-efficient scheme. Using exact exchange and the RPA 
correlation energy yields improved bulk moduli and lattice constants (e.g. a mean absolute 
error in lattice constants of 0.6% compared to 1.2% for the LDA). Apart from Na, the 
calculated RPA lattice constants are larger than the experimental values, a result also found 
in Ref. 46. 
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TABLE I. Lattice constants (in A) and bulk moduli (GPa) calculated for the test set of 10 materials compared to the 
experimental data tabulated in Ref. 46. The results were calculated from correlation energies obtained by the extrapolation 
procedure of equation 39. The mean absolute error (M.A.E.) compared to experiment is shown in the final row. The experi¬ 
mental lattice constants were corrected for expansion due to zero-point motion; the bulk moduli have not been corrected. For 
comparison the LDA and RPA calculations of Ref. 46 are also presented (note that these RPA calculations were performed on 
top of Kohn-Sham states obtained within the generalized-gradient approximation^^). The CP and JGMs kernels coincide for 
metallic systems. 


LDA LDA^ PBE RPA RPA*^ rALDA rALDAc CDOP CDOPs CP CPd JGMs Exp. 


c 

3.533 

3.534 

3.573 

3.566 

3.572 

3.563 

3.562 

3.561 

3.565 

3.562 

3.553 

3.550 

3.553 

Si 

5.407 

5.404 

5.477 

5.449 

5.432 

5.456 

5.453 

5.446 

5.452 

5.454 

5.450 

5.437 

5.421 

SiC 

4.338 

4.332 

4.390 

4.380 

4.365 

4.380 

4.379 

4.374 

4.379 

4.378 

4.371 

4.361 

4.346 

MgO 

4.163 

4.169 

4.255 

4.229 

4.225 

4.229 

4.228 

4.224 

4.231 

4.228 

4.222 

4.202 

4.189 

LiCl 

4.972 

4.967 

5.157 

5.082 

5.074 

5.099 

5.097 

5.069 

5.094 

5.093 

5.095 

4.996 

5.070 

LiE 

3.926 

3.913 

4.080 

4.010 

3.998 

4.011 

4.011 

3.989 

4.013 

4.010 

4.007 

3.978 

3.972 

A1 

3.987 

3.983 

4.044 

4.042 

4.037 

4.053 

4.051 

4.041 

4.048 

4.050 

4.051 

4.050 

4.018 

Na 

4.054 

4.056 

4.197 

4.205 

4.182 

4.229 

4.225 

4.221 

4.233 

4.232 

4.257 

4.232 

4.214 

Cu 

3.530 

3.523 

3.643 

3.622 

3.597 

3.612 

3.616 

3.622 

3.625 

3.625 

3.626 

3.625 

3.595 

Pd 

3.839 

3.830 

3.941 

3.914 

3.896 

3.919 

3.918 

3.916 

3.920 

3.920 

3.918 

3.920 

3.876 

% M.A.E. 

1.2 

1.3 

1.3 

0.6 

0.5 

0.7 

0.7 

0.5 

0.7 

0.7 

0.7 

0.6 

_ 


LDA LDA'^ PBE RPA RPA'^ rALDA rALDAc CDOP CDOPs CP CPd JGMs Exp. 


C 

465 

465 

432 

435 

441 

439 

440 

443 

439 

439 

448 

455 

443 

Si 

95 

97 

88 

95 

99 

95 

95 

96 

95 

95 

99 

98 

99 

SiC 

228 

229 

211 

220 

223 

220 

221 

223 

221 

221 

224 

231 

225 

MgO 

172 

172 

151 

163 

168 

161 

161 

164 

161 

162 

163 

174 

165 

LiCl 

41 

41 

31 

36 

37 

35 

35 

38 

36 

36 

35 

40 

35 

LiE 

86 

87 

67 

71 

76 

71 

71 

77 

72 

72 

72 

72 

70 

A1 

83 

84 

77 

78 

77 

76 

76 

78 

76 

76 

75 

76 

79 

Na 

9 

9 

8 

8 

8 

8 

8 

8 

8 

8 

7 

8 

8 

Cu 

184 

186 

136 

156 

153 

156 

156 

155 

151 

151 

152 

151 

142 

Pd 

226 

226 

169 

199 

181 

195 

196 

196 

194 

194 

192 

194 

195 

% M.A.E. 

12 

12 

7 

3 

3 

3 

3 

4 

3 

3 

3 

4 

_ 


Ref. 46 


25 















3. rALDA and rALDAc 


For the non-metallic systems, the rALDA and rALDAc kernels prodnce lattice constants 
and bnlk modnli which are essentially indistinguishable from each other. In turn, these 
results are in close agreement with the RPA. For metals, one can identify differences between 
the kernels, although the magnitude of variation is still very small (<0.1%). The close 
agreement between rALDA and rALDAc conhrms that the exchange contribution dominates 
in equation 10, and supports the use of the exchange-only rALDA kernel. 

The rALDA kernels also display the fastest convergence with respect to i^cut- Recalling 
the form of the kernels (equation 14), the components of Jhxc are truncated for wavevectors 
exceeding the cutoff kc- For a homogeneous system (equation 8), for k' > kc, the interact¬ 
ing and non-interacting response functions coincide and therefore the contribution to the 
correlation energy at these wavevectors vanishes. In inhomogeneous systems, high-density 
regions (large kc) give terms that converge like the RPA, but the rALDA convergence is still 
superior after the kernel averaging procedure (equation 38) is applied. 

From this convergence behavior, we conclude that the short-range description of correla¬ 
tion obtained from an XC-kernel like the rALDA is easier to describe in a plane-wave basis 
than the erroneous short-range behavior of the RPA. This result might have been anticipated 
from the HEG, where the coupling-constant averaged pair-distribution function calculated 
for the rALDA is softer than for the RPA (Fig. 2). 


4- CP and CDOPs 

The CP and CDOPs kernels yield lattice constants and bulk moduli which are also very 
similar to each other across the full range of systems. This behavior can be explained by 
considering Fig. 1(a), where it can be seen that lies more negative than for k less 

that ~ ‘2kF, and more positive otherwise. The kernel averaging procedure smears out these 
differences. In particular, there is a negligible effect from modifying the large-A; behavior 
from —4:7rB/k‘^ (CDOPs) to —47r//c^ (CP). 

and closely follow the rALDA kernels (and the RPA) for non-metallic systems. 
For metallic systems, differences of order 0.3% can be observed. The most likely reason for 
this difference is the long-range behavior of the rALDA kernels, which display decaying 
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oscillations, compared to the CDOPs and CP kernels which go to zero more smoothly. The 
small positive hump displayed by the CDOPs kernel in real space [Fig. 1(b)] appears to have 
little effect on the correlation energy. 


5 . CD OP 

The CDOP kernel (equation 17) differs from by having a local term. This local 

term has a noticeable effect on the calculated structural properties, with the CDOP kernel 
predicting slightly smaller and larger lattice constants and bulk moduli, respectively. Indeed 
the CDOP kernel displays the closest agreement with experimental lattice constants, but 
performs less well on bulk moduli. 

In section IIC 4 it was pointed out that the local term in the CDOP kernel leads to 
a divergent pair-distribution function. The local term may also be expected to introduce 
convergence problems, as demonstrated for the (entirely local) ALDA kernel.In the 
current work, we have not found any signihcant difference in the convergence behavior of 
the CDOP and CDOPs kernels when calculating lattice constants and bulk moduli for Ecut < 
400 eV. Only in cases where the RPA correlation energy converges relatively quickly (e.g. Al) 
can we observe a slowly-converging positive contribution to the CDOP correlation energy 
which is reminiscent of that found for the ALDA, c.f. Fig. 3 of Ref. 30. However, unlike 
for the ALDA, the magnitude of this contribution is very small compared to the RPA-like 
convergence (e.g. Fig. 3). 


6. CPd 

The dynamical CPd kernel displays slight differences to its ca = 0 limit, Compared 
to the static kernels where the range of f^c is fixed by the density, for the CPd kernel the 
frequency appearing in the denominator of xks also affects the range. Interestingly the CPd 
bulk moduli of insulators are slightly closer to experiment. In other cases we hnd that the 
CP and CPd kernels predict similar results except for Na, where the CPd kernel hnds a 
larger lattice constant and smaller bulk modulus, and C, where the CPd lattice constant lies 
on top of the experimental value. 

The CPd results show that even a simple dynamical kernel can predict different struc- 
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TABLE II. Parameters relating to the JGMs kernel. The values of {a) were obtained by inserting 
the experimental band gaps Eg into equation 40, while inserting the “effective gaps” E®® yields 
the aLRc values reported in Ref 70 (these calculations were performed at the experimental lattice 
constant). The olrc values for LiCl and LiF were obtained from equation (4) of Ref. 70 using the 
dielectric constants tabulated in Ref 91. The experimental (direct) band gaps were obtained from 
Refs. 92-97. 



c 

Si 

SiC 

MgO 

LiCl 

LiF 

Eg (eV) 

7.3 

3.4 

6.0 

7.8 

9.4 

14.2 

(a) 

0.58 

: 0.89 

1 1.30 

2.32 

5.60 

7.03 

Ef (eV) 

7.43 

. 1.57 

■ 3.62 

6.74 

3.98 

6.07 

OLRC 

0.6 

0.2 

0.5 

1.8 

1.5 

2.2 


tural properties. This result is not obvious from studies on the HEG, where tests on the 
more complicated frequency-dependent kernel of Ref. 18 found dynamical effects to be less 
important than nonlocality when calculating correlation energies. Although we do not ob¬ 
serve systematic improvement with the CPd kernel, it would be interesting to investigate its 
performance for systems with a greater degree of inhomogeneity, e.g. molecules and surfaces. 


7. JGMs 


The lattice constants calculated with the JGMs kernel for insulating systems display the 
closest agreement with experiment out of all of the considered kernels, except for the notable 
example of LiGl, where the JGMs lattice constant is underestimated by 1.4%. However the 
agreement with experimental bulk moduli is poorer, in some cases (SiG, MgO) worse than 
the LDA. For metallic systems, the JGMs and GP kernels coincide. 

It is important to establish the importance of the value of Eg. In the current work we have 
used the experimental, direct gap (Table II), but we equally could have chosen the indirect 
gap, or even dehned a more general r-dependent gap function.®® An alternative option is 
to make the link to the description of excitons®®’^® and consider the head (G = G' = 0) of 
in the g —>■ 0 limit, which can be written as —{a)/q^ where 


(«) 


dr 
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FIG. 6. Percentage deviations of lattice constants compared to experiment,^® calculated with the 
RPA, and the JGMs kernel using the experimental direct band gaps or effective LRC gaps (see 
Table II). 

The values of (a) computed from equation 40 for the experimental gaps are given in Table II. 
These values can be compared to Ref. 70, where a long-range (LRC) attractive kernel was 
introduced as fxc{R) = — Q;LRc/(4vri?). We note that the head of this matrix in reciprocal 
space in the g —>■ 0 limit coincides with the JGMs kernel with olrc and also that 

this single matrix element is considered the most important for the calculation of excitonic 
effects.™ 

From Table II it is clear that the values of (a) calculated with the experimental gaps and 
PAW densities are somewhat larger than the values of olrc reported in Ref. 70, which were 
found to give a good description of excitonic effects in absorption spectra of semiconductors 
and MgO. To explore this point further we adopted an inverse approach and considered an 
effective gap Eg^, which when inserted into equation 40 yields arRC- These LRC “gaps” are 
smaller than experimental values, especially for the ionic compounds. Indeed the empirical 
Q^lrc values of LiCl and LiF are signihcantly smaller than those expected both from the 
JGM or bootstrap kernels,which have been shown to accurately capture the exciton in 
LiF. 

We repeated the JGMs kernel calculations using the LRC gaps Eg^, and show the obtained 
lattice constants in Fig. 6. The LRC results he between the lattice constants calculated with 
the RPA and with the JGMs kernel/experimental gaps, and thus improve the LiCl result. 
Comparison of LiCl and LiF demonstrates the nonlinear relation between Eg^ and the lattice 
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constant. In both cases the effective gap is reduced by more than 50% from its experimental 
value, but the effect on the LiCl lattice constant is an order of magnitude larger than for 
LiF. 

The improved agreement of lattice constants with experiment compared to the RPA shows 
that XC-kernels with long-range components represent an interesting avenue to study. A key 
question is whether the tendency for the JGMs kernel to favor smaller lattice constants than 
the RPA is directly related to the fundamental long-range character of the former, or whether 
it is in fact a consequence of the precise form of the kernel. The strength of the long-range 
part of the JGMs Hartree-XG kernel is determined by exp(—A'^/47rn) (equation 23), which 
becomes RPA-like in the high density limit and vanishes in the low density limit. The RPA 
correlation energy is generally negative, while a zero Jhxc implies a zero correlation energy. 
Interpolating these two limits implies that a more negative (i.e. stable) JGMs correlation 
energy will correspond to a higher density, thus favoring a lower lattice constant. This 
observation also provides an explanation for the varying behavior of the bulk modulus and 
also the strong nonlinearity in the variation of the lattice constant with band gap, since the 
energy-volume relation is expected to be sensitive to the relative magnitude of Eg and n. 

We note that the bootstrap approach^^ is an alternative method of constructing a long- 
range kernel. Since the bootstrap kernel is constructed from xkS; using it would avoid both 
the input of Eg and the averaging procedures discussed in section IIG. However, it would 
be necessary to ensure that the bootstrap kernel displayed reasonable behavior in the large 
/c-limit. 

B. Absolute correlation energies 

In Fig. 7 we show the absolute correlation energy per electron calculated using each 
of the different kernels for the materials in the test set. Absolute correlation energies are 
generally considered less robust than properties constructed from energy differences, being 
more difficult to converge and sensitive to details of the PAW potentials. However one can 
still perform a comparison between kernels, and look for similarities with the trends observed 
for the HEG [Fig. 1(d)]. 

The most obvious feature of Fig. 7 is the reduction of absolute correlation energy on 
moving from the RPA to a nonzero fxc, ranging from 0.1 eV for Na to 0.5 eV for Si. This 
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FIG. 7. Absolute correlation energies Ec (equation 1) calculated per (valence) electron for different 
kernels. The correlation energy obtained for Si from the diffusion Monte Carlo (DMC) calculations 
of Ref. 99 is also shown. The core/valence partitioning of the PAW potentials is given in section IIH. 
The calculations were performed at the experimental lattice constant. 


change is the same order of magnitude as observed for the HEG. The ordering of the HEG 
correlation energy calculated with different kernels is also largely preserved, with the GP 
and GDOPs kernels predicting larger magnitudes than GDOP and the rALDA kernels. 

The difference between the rALDA and rALDAc kernels is small, with the rALDA cor¬ 
relation energy being more negative by order 1% or 0.01 eV per electron. The difference 
between the static and dynamical forms of the GP kernel is an order of magnitude larger, 
with the static correlation energy being more negative. Meanwhile the removal of the local 
term in the GDOP kernel increases the magnitude of the correlation energy, with the GDOPs 
having a more negative correlation energy than GDOP by 5% or 0.05 eV per electron. 

As in Ref. 29 we can tentatively compare our calculated correlation energy for Si with 
the diffusion Monte Garlo (DMG) calculations of Ref. 99. Reassuringly the DMG correlation 
energy lies among the values calculated with the model exchange kernels (Fig. 7), in fact 
displaying closest agreement with GPd, rALDA and GDOP kernels. We also note that our 
calculated GDOP correlation energy for Si (-1.02 eV per electron) lies on top of the value 
recently reported in Ref. 29 using a pseudopotential approximation and a different averaging 
scheme (equation 37). With highly-accurate calculations of correlation energies in extended 
systems now becoming a reality, comparisons of this sort should become a useful test for 
new kernels. 
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C. Kernel averaging scheme 


It is interesting to compare the structural properties and correlation energies calcu¬ 
lated using the symmetrized-wavevector averaging scheme (equation 38) to the two-point 
symmetrized density (equation 34). The latter has previously been implemented for the 
rALDA,^^ so here we restrict the comparison to this kernel. 

Considering the lattice constants hrst, we typically hnd a difference of 0.2% between the 
two methods, with the symmetrized-density values larger than those calculated with the 
symmetrized-wavevector in most cases. Interestingly the agreement is worse for the bulk 
moduli, with an average deviation of 6%. The absolute correlation energies also show a 
larger (4%) deviation, where using the two-point symmetrized density scheme consistently 
yields more negative rALDA correlation energies than the symmetrized-wavevector scheme, 
by an average of 0.04 eV per electron. 

To understand the origin of these differences it is necessary to consider the practical 
implementation of the two-point density average (equation 34). As mentioned in section IIG, 
constructing the kernel in this way involves sampling the 1/R Coulomb interaction in real 
space. The divergence at i? = 0 is replaced with a spherical average oi 1/R taken over 
the volume per point in the real-space grid used to evaluate the integral.The absolute 
correlation energy is therefore rather sensitive to this grid spacing, and its dependence on 
volume (i.e. the bulk modulus) will also be difficult to converge. 

The symmetrized-wavevector approach only samples the density on the real space grid, 
and therefore shows a much weaker dependence on the spacing between the grid points. We 
verified this behavior for diamond C, where the symmetrized-wavevector correlation energy 
changes by less than 10“^ eV/electron on reducing the grid spacing from 0.17 to 0.11 A. This 
is several orders of magnitude faster than the symmetrized-density approach,illustrating 
a computational advantage of equation 38. 


D. Spin and atomization energies: the H 2 molecule 

Throughout this study we have not considered any spin-dependence of the XC-kernels. 
However, the calculation of atomization or cohesive energies usually requires the description 
of spin-polarized atoms or molecules. In this section we provide a demonstration of the 
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TABLE III. Correlation energies of H, H 2 and He and the atomization energy of the H 2 molecule 
[£'at(H 2 )] calculated at different levels of theory. The CCSD values are taken from Ref. 101 and 
the experimental atomization energy from Ref. 102. The rALDA results were calculated including 
spin-polarization, with the kernels constructed either from equation 38 or equation 34. All values 
are given in eV. 


RPA rALDA^^ rALDAi^ CCSD Exp 


H 

-0.57 

0.06 

-0.02 

0.00 — 

H 2 

-2.22 

-1.04 

-1.22 

-1.11 — 

He 

-1.82 

-1.00 

-1.08 

-1.14 — 

Aat(H2) 

4.74 

4.82 

4.85 

4.75 4.75 


^ Symmetrized wavevector, equation 38 
Symmetrized density, equation 34 

importance of spin by calculating the atomization energy of the H 2 molecule with the rALDA 
kernel. 

First, we note that the symmetrized-wavevector averaging procedure (equation 38) can 
be equally applied to extended and hnite systems. In fact, for the rALDA kernel one 
can exploit the fact that the Hartree-XC kernel strictly vanishes at any points in space 
where the density is less than (|G -|- q||G' -|- q|)^/^/(247r^). Therefore the Fourier transform 
can be performed in a small box which excludes the vacuum region generally required to 
model isolated systems with periodic boundary conditions. Since the H 2 molecule is spin- 
unpolarized we can calculate its correlation energy without any further consideration, and 
obtain a value of -1.04 eV with the rALDA kernel. This value is within 0.1 eV of the value 
obtained from coupled-cluster calculations^®^ and a signihcant improvement (>1 eV) over 
the RPA. We hnd a similar level of agreement for the He atom (Table III). 

For the spin-polarized H atom, following the analysis of Ref. 31 we replace the integral 
equation (equation 3) with its the spin-polarized version, valid for systems where only one 
spin channel is occupied: 

x^^(q, uj) = A]<s(q) + XKs(q> ^)/lrL(q, w)A^^(q, (4i) 

The above quantities are related to equation 1 through the simple substitutions y —)■ 
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and xks Xks- 

Proceeding fnrther reqnires the spin-polarized form of the Hartree-XC kernel 
To onr knowledge, of the XC-kernels stndied in this work has been derived only for the 
rALDA, given as:^^ 



2 X e{kc 



+ 6{k — kc) 


dTT 


(42) 


Eqnation 42 differs from the spin-nnpolarized rALDA expression (eqnation 14) by a factor 
of two in front of the part of the kernel corresponding to the ALDA, reflecting the fact 
that the exchange interaction acts only between electrons with the same spin. Using eqna- 
tions 41 and 42 to calculate the correlation energy of the H atom yields a value of 0.06 eV, 
compared to the exact value of zero and an RPA value of -0.57 eV. 

Taking the H and H 2 calculations together yields an rALDA atomization energy of 
4.82 eV, which is within 0.1 eV of the experimental value of 4.75 eV.^°^ We note that 
the RPA benefits from substantial error cancellation and yields an atomization energy very 
close to experiment (4.74 eV, Table III). However the H 2 molecule is a rather special case, 
and the RPA usually demonstrates percentage errors of order 10% in atomization energies. 
The rALDA kernel corrects the correlation energies of the individual H 2 and H systems and 
maintains close agreement with the experimental atomization energy. 

In Table III we also present the rALDA correlation energies using the two-point density 
average, equation 34. As found in bulk systems, the correlation energies calculated with the 
two-point density average are more negative (~ 0.06 eV/electron) than those calculated with 
the symmetrized wavevector. However the agreement in atomization energies is better than 
0.03 eV. We hnd it encouraging that the symmetrized-wavevector approach gives such similar 
results to the more intuitive two-point density average when calculating the atomization 
energy. 

We note that if we do not use the spin-polarized form of the kernel (equation 42), we 
hnd a correlation energy of -0.17 eV for the H atom and an atomization energy of 4.37 eV. 
This value is in signihcantly worse agreement with experiment than the RPA or even the 
LDA (4.89 eV), emphasizing the importance of a rigorous treatment of spin. An important 
direction for further study is the introduction of spin-dependence into kernels derived from 
the spin-nnpolarized HEG. 
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IV. CONCLUSIONS 


We have calculated the correlation energy of a test set of 10 materials within the adiabatic- 
connection fluctuation-dissipation formulation of density-functional theory (ACFD-DFT). 
We used a hierarchy of approximations for the exchange-correlation (XC) kernel fxc, includ¬ 
ing the random phase approximation (RPA, f^c = 0), the recently-introduced renormal¬ 
ized kernels (rALDA),^^ a kernel which satishes the exact static limits of the electron gas 
(CDOP)/^ a model dynamical kernel (CPd)^® and a kernel which diverges oc in the 
smalh/c limit (JGMs).^® In order to apply homogeneous kernels to inhomogeneous systems 
we applied a reciprocal space averaging scheme employing wavevector symmetrization.^® 
For each kernel and material pair we calculated the lattice constant and bulk modulus, and 
compared our results to previous calculations and experiment.^® 

For all materials, including a nonzero fxc reduces the magnitude of the correlation energy 
compared to the RPA by 0.1-0.5 eV per electron. This result mirrors the homogeneous 
electron gas (HEG), where the RPA correlation energy is too negative by at least 0.3 eV 
over a wide range of densities.^® However the variation in correlation energy between each 
fxc is much smaller, on the scale of 0.01-0.1 eV per electron. Encouragingly the correlation 
energies calculated with XG-kernels are found to he very close to diffusion Monte Garlo 
data available for Si.®® Furthermore, calculations with XG-kernels display faster basis-set 
convergence than the RPA due to the suppression of high energy plane-wave components of 
the Goulomb potential. 

Gonsidering lattice constants and bulk moduli, we found only small variations between 
the RPA and different XG-kernels. In particular, static XG-kernels that only satisfy the 
/c —)■ 0, a; = 0 limit of the HEG (rALDA, GP, GDOPs) produce very similar results. The 
structural properties calculated with the dynamical GPd kernel are in better agreement with 
experiment in some cases (e.g. the bulk moduli of non-metallic systems), but the improve¬ 
ment is not systematic (e.g. Na). Satisfying the k —)■ oo, u = 0 limit of the HEG (which 
adds a local term to fxc, e.g. the GDOP kernel) also yields good agreement with experi¬ 
mental lattice constants, despite the kernel having a diverging pair-distribution function.^® 
Finally, the JGMs kernel predicts a reduction in lattice constants and an increase in bulk 
moduli for non-metallic systems, bringing the former into closer agreement with experiment. 
The current study however cannot distinguish whether this behavior is due to the general 
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long-range —Q;/(47ri?) character of the kernel/'^ or to the density-dependence of a specihc 
to the JGMs model. 

The ACFD-DFT scheme described here clearly involves a number of choices, including (a) 
the approximation used to generate the noninteracting response function XkS) (b) the form 
of fxc (including spin-dependence) and (c) the averaging scheme used to generalize a HEG 
XG-kernel to an inhomogeneous system. Fixing all factors except (b), as we have done here, 
points us towards the essential properties of a model fxc- Based on the similar performance of 
the different XG-kernels, the current work supports the idea that fxc should be kept as simple 
as possible, i.e. be static, tend to a density-dependent constant at small k and decay oc 
at large k. In this respect the exchange-only rALDA kernel is attractive, since it scales 
simply with the coupling constant A and has good convergence properties. The introduction 
of additional computational expense and uncertainty associated with a dynamical kernel, a 
divergence oc l//c^ at /c = 0 or even a local term in fxc is difficult to justify based on the 
performance of the GPd, JGMs and GDOP kernels for lattice constants and bulk moduli, 
although each kernel was found to offer improved agreement with experiment in certain 
cases. 

On the other hand, by focusing on the structural properties of bulk solids we have cho¬ 
sen systems where the RPA already performs very well. It is encouraging that the model 
XG-kernels can maintain this good performance whilst correcting the magnitude of the cor¬ 
relation energy by several eV per atom, but arguably their real test lies in cases where the 
RPA is less successful. Already the rALDA kernel has been found to improve the description 
of atomization and cohesive energies^^’^^ but a number of challenges remain, particularly in 
the description of molecular dissociation.The framework described in the current 
study provides the base for the application of a full range of kernels to these more challenging 
systems. 
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